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Abstract 


Two  vertical  soft  landing  problems  are  investigated  in  this  report.  First,  the  soft  landing 
problem  of  quadrotor  in  presence  of  unmodelled  dynamics  is  investigated.  A  neural  network 
based  disturbance  estimation  is  adopted  to  capture  the  unmodeled  quadrotor  dynamics  due  to 
rotor  blade  flapping  phenomenon.  An  adaptive  guidance  law  with  the  Dynamic  Inversion  (DI) 
as  baseline  algorithm  is  illustrated  for  soft  vertical  touch  down.  Next,  the  autonomous  landing 
of  a  spacecraft  on  the  lunar  surface  is  explored.  To  ensure  the  smooth  touchdown  of  the 
spacecraft  on  the  lunar  surface,  a  nonlinear  optimal  control  theory  based  Generalized  model 
predictive  static  programming  (G-MPSP)  guidance  is  proposed.  As  the  G-MPSP  formulation 
incorporates  the  terminal  condition  as  a  hard  constraint,  it  ensures  the  high  terminal  accuracy 
of  position  and  velocity  of  the  spacecraft.  Also  the  vertical  orientation  of  the  spacecraft  during 
touchdown  is  achieved  through  the  soft  constraint  formulation  by  the  proper  selection  of  the 
control  weight  matrix.  Effectiveness  of  the  proposed  guidance  methods  are  demonstrated 
using  simulation  results. 
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0.1  Introduction 

In  general  the  objective  of  the  guidance  algorithm  for  the  autonomous  vertical  soft  landing  of 
a  UAV  is  to  ensure  that  it  touches  down  at  the  designated  site  vertically  with  near-zero  vertical 
velocity.  Near-zero  touchdown  velocity  is  essential  for  the  safety  of  the  vehicle  as  well  the 
on-board  payload/instrument.  The  desired  lading  site  comes  from  the  mission  requirement. 
During  landing  it  is  also  essential  to  maintain  the  vertical  orientation  of  the  vehicle  with  re¬ 
spect  to  ground.  To  avoid  the  local  uneven  surface  of  desired  landing  site  it  might  be  necessary 
to  shift  the  landing  site  location  to  a  nearby  smooth  ground.  Such  requirements  may  comes 
from  the  vision  based  information  and  the  necessary  correction  for  the  trajectory  needs  to  be 
recomputed  online.  Two  related  but  different  problems  (which  offer  different  challenges)  are 
investigated  in  this  report.  One  problem  deals  with  the  quadrotor  UAV  which  needs  to  land  on 
earth  surface  in  presence  of  modelling  uncertainty  and  another  deals  with  soft  landing  on  the 
surface  of  the  moon. 


Quadrotor  UAVs  are  emerging  as  popular  choice  as  both  experimental  research  platform 
and  for  commercial  use.  Due  to  their  design  and  construction  simplicity,  vertical  take  off 
and  landing,  hovering  capability  and  ease  of  operation  within  confined  spaces,  capability  to 
do  aggressive  maneuvers  quadrotors  are  preferred  over  winged  UAV  Quadrotor  utilizes  lift 
generated  from  four  propellers  blades  for  the  lateral,  longitudinal  and  altitude  control.  Unlike 
conventional  helicopters  which  utilizes  a  tail  rotor  for  controlling  body  spin  resulting  from 
angular  momentum  conservation  of  the  rotor  blades,  quadrotor  does  not  need  a  tail  rotor 
since  opposite  rotor  blade  pair  rotate  in  opposite  direction  and  hence  counteracting  each 
others  angular  momentum,  this  feature  makes  quadrotors  yet  more  reliable  over  helicopters. 
For  detailed  operation,  functional  details  of  quadrotor  system  one  can  refer  Beard  [2008]. 


After  1963,  the  first  manned  flight  (hovering)  of  the  four  rotor  configuration  of  quadrotor 
helicopter  by  Cutis-Wright  X19A  no  much  advancement  was  seen  in  the  multi-rotor  flight 
due  to  stability  issues  due  actuator  dissimilarities  which  lead  to  instability  in  the  quad  flight, 
rendering  stable  flight  nearly  impossible.  With  recent  advancement  of  the  micro-controller 
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and  solid  state  electronics  the  onboard  computer  has  become  more  smaller  and  faster  Hart¬ 
ley  et  al.  [2013].  With  availability  of  small,  light,  high  fidelity  sensors  (Inertial  Measurement 
Units  IMU)  and  processors  on  board,  advanced  controls  for  quadrotor  have  become  a  reality. 
Many  research  teams  have  proposed  various  control  architectures  for  quadrotor  control  and 
guidance.  Most  popular  control  experimented  on  quadrotor  UAV  platform  are  linear  feedback 
control  architecture  like  PID,  LQR  Huang  et  al.  [2009],  Pounds  et  al.  [2010]  etc.  For  indoor 
and  low  speed  flight  fixed  gain  linear  feedback  controller  are  quite  capable  of  delivering  the 
requisite  performance  of  attitude  control  and  trajectory  tracking.  Since  linear  control  utilizes 
the  linearized  plant  dynamics  model  for  control  computation,  under  scenario  like  parameter 
uncertainty  viz  inertia  change  or  partial  actuator  failure,  under  influence  of  un-modeled  plant 
dynamics  or  under  external  perturbation  these  controller  fail  to  achieve  the  required  perfor¬ 
mance.  Many  research  group  have  demonstrated  the  command  trajectory  tracking  for  indoor 
flights  and  formation  flight  like  OS4  quadrotor  project  Bouabdallah  et  al.  [2005],  MIT  SWARM 
Valenti  et  al.  [2006]  and  at  Stanford  Testbed  of  Autonomous  Rotorcraft  for  Multi -Agent  Control 
(STARMAC)  Huang  et  al.  [2009] ,  Hoffmann  et  al.  [2007] .  Several  adaptive  control  techniques  are 
also  experimented  as  well  in  quadrotor  control,  Annaswamy  et  al.  have  experimented  MRAC 
based  controller  in  presence  of  actuator  uncertainties  Dydek  et  al.  [2010]  and  Whitehead  and 
Bieniaswski  Whitehead  and  Bieniawski  [2010]  have  demonstrated  MRAC  controller  for  step 
command  in  altitude  tracking  under  actuator  degradation.  Chowdhary  et  al.  Chowdhary  et  al. 
[2012]  have  experimented  and  presented  concurrent  learning  MRAC  for  outer-loop  control 
and  PD  based  controller  for  attitude  controller  of  quad-rotor  where  the  gains  are  tuned  for 
larger  quad-rotor  frame  and  same  control  gains  are  ported  to  much  smaller  quad-rotor  for 
trajectory  tracking  problem.  Gabriel  Hoffman  et  al.  have  experimented  PD  controller  with 
integral  feedback  controller  catering  for  blade  flapping  uncertainty.  Integral  feedback  control 
accounts  for  this  uncertainty  as  a  constant  bias  i.e.  for  uncertainty  at  constant  velocity  and 
attitude  angle,  but  as  the  vehicle  velocity  and  attitude  changes  the  integral  control  has  to  adapt 
to  new  uncertain  moment  values.  In  the  problem  experimented  in  this  paper  for  commanded 
spiralling  trajectory  tracking  this  control  architecture  will  fail  to  perform  the  desired  tracking, 
hence  there  is  need  for  controller  which  adapts  online  to  approximate  the  uncertainty  for  con¬ 
trol  computation.  Similarly  Pounds  et  al.  Pounds  et  al.  [2010]  have  experimented  linear  PID 
and  LQR  based  controller  compensating  for  uncertainty  due  to  blade  flapping  phenomenon. 
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But  as  per  authors  knowledge  not  many  have  experimented  on  adaptive  control  catering  for 
aerodynamic  uncertainty  in  trajectory  tracking. 


This  report  proposes  a  control  architecture  utilizing  the  well-established  Dynamic  Inver¬ 
sion  (DI)  theory  and  augmenting  it  with  online  trained  Radial  Basis  Function  (RBF)  neural 
networks  a  robust  nonlinear  controller  catering  to  the  actual  plant  model  is  experimented  and 
results  are  presented  in  this  article.  This  control  architecture  is  experimented  with  plant  model 
of  Quadrotor  trajectory  tracking  problem.The  RBF-NN  is  used  to  capture  the  unmodeled 
dynamics /uncertainty  in  the  system  due  to  rotor  blade  flapping  phenomenon.  A  new  weight 
update  rule  based  adaptive  controller  with  DI  as  baseline  controller  is  validated  for  quadrotor 
trajectory  tracking  problem.  It  is  observed  that  proposed  controller  is  able  to  capture  the 
time  varying  uncertainty  perfectly  and  thereby  augmenting  the  approximate  system  dynamics 
to  track  the  actual  system  dynamics.  The  DI  control  evaluated  using  online  augmented  ap¬ 
proximate  plant,  satisfactorily  tracks  the  commanded  attitude  angles  and  hence  achieving 
the  desired  trajectory  tracking  through  outer-loop  control.  The  details  are  presented  in  the 
subsequent  results  section. 
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Control  Synthesis:  Generic  Theory 


In  this  section  the  methodology  used  for  autonomous  guidance  for  soft  landing  is  described 
with  necessary  mathematical  details. 


1 . 1  Philosophy  of  Adaptive  Dynamic  Inversion 

This  section  gives  the  details  of  the  philosophy  of  the  online  adaption  method  of  DI  controller. 
A  online  trained  NN  is  used  along  with  the  inner  loop  DI  controller  of  the  quad-rotor  to  capture 
the  plant  disturbance,  update  the  plant  model  and  there  by  synthesizing  the  over  all  controller. 
RBF-NN  approximates  the  un-modeled  dynamics  which  is  crucial  for  online  adaptation  of 
DI  controller.  The  approximate  model  synthesized  using  trained  NN,  is  termed  as  a  modified 
state  observer  Joshi  and  Padhi  [2013].  NNs  use  the  channel  wise  error  information  in  the  state 
for  training,  single  layer  networks  are  employed  with  radial  basis  function  as  basis  vectors. 
The  nominal  plant  model  which  is  fully  known  is  assumed  to  have  the  following  structure 

X  =  foiX)  +  g(X)U  (1.1) 

where  X  e  Jf";  U  e  Jfra  and  y  e  5fp,  fo(X)  and  g(X)  represents  nominal  system.  Nominal 
controller  U  is  designed  such  that  the  states  of  nominal  plant  model  X  track  the  respective 
states  of  the  desired  plant  (1.4)  i.e.  goal  is  to  ensure  that  X  —  Xj  as  t  —  oo. 
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The  actual  system  or  actual  plant  model  is  considered  as  to  be  having  a  structure 

X  =  f(X)  +  g(X)U+dext(X)  (1.2) 

where  f(_X)  is  actual  system  model,  dext(X)  denotes  the  disturbance  external  to  the  system. 
Since  the  actual  system  parameters  are  unknown  and  external  disturbance  model  is  uncertain, 
the  actual  plant  model  cannot  be  used  for  control  synthesis. 

The  total  system  uncertainty  can  therefore  be  denoted  by  the  term  d{X)  =  ( f(X)  -  fo(X)  + 
dext(X)),  where  d(X)  e  -Rn  represents  the  total  uncertainty  term  in  the  system  (1.2).  The 
uncertain  actual  plant  model  can  be  represented  using  nominal  plant  model  and  uncertainty 
term  d{X)  as  follows 


X  =  f0{X)  +  giX)U  +  d[Xy,  y=h(X)  (1.3) 

Consider  a  nonlinear  systems  representation  for  the  desired  plant  dynamics.  The  desired 
plant  dynamics  output  the  desired  state  trajectory  for  the  actual  system  to  follow 


Xd  =  fd(Xd) 


(1.4) 


where  Xd  e  5R"  represents  the  desire  state  vector. 

In  the  actual  plant  model  (1.3)  the  term  d[X)  is  unknown,  the  objective  is  to  first  capture  the 
unknown  function  d{X)  through  NN  approximation  d{X)  =  14/7  0(X)  and  form  an  approxi¬ 
mate  system  model  as  follows. 

=  f0 (X)  +  g(X)U+  d(X)  +  KAX-  Xa)  (1.5) 

Concept  of  modified  state  observer  is  applied  with  DI  to  form  the  actual  controller.  The 
controller  synthesis  call  for  two  step  process  for  control  realization 

1.  X  —  Xa  as  t  —  oo.  The  complete  state  feedback  is  not  available  but  reduced  order  state 
vector  Xmeas  is  available  through  the  state  measurement.  The  measured  state  values 
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are  considered  to  be  corrupted  by  measurement  noise  and  can  be  represented  as 

Xmeas  =  T  +  (0,  O'  ) 

where  o  represents  the  standard  deviation  of  the  zero  mean  white  noise.  Hence  the 
revised  the  aim  is  to  ensure 

Xest  Xa 

where  Xest  is  the  Extended  Kalman  filter  estimate  of  the  state  vector  X.  Using  estimated 
state  value,  the  approximate  plant  model  is  written  as 


Xa  =  MXest)  +  g(Xest)U+diXest)  +  KT{Xest  -  Xa)  (1.6) 

2.  Xa  —  Xa  as  t  — •  oo:  This  process  is  accomplished  through  control  synthesis  using  online 
adaptation  of  DI  controller. 

1.1.1  Control  Synthesis  using  Feed-Back  Linearization  Method 

For  the  control  synthesis  it  is  considered  that  the  unknown  function  d(Xest)  is  function  of  Xest 
and  not  X  since  the  uncertainty  approximation  available  through  RBF-NN  and  the  RBF  basis 
is  function  of  Xest. 

®(Xesf)  =  e^Xest~C^  (1.7) 

With  the  objective  of  the  controller  to  enforce  Xa  -*■  X,i.  The  plant  dynamics  for  Feedback 
Linearization  control  synthesis  can  be  written  as  follows, 

Xa  =  v  (1.8) 

where  v  is  termed  as  synthetic  control  to  linearized  plant  model  (1.8).  Any  linear  control 
technique  can  be  implemented  to  evaluate  the  linear  synthetic  control  term  v.  In  this  paper 
a  Linear  Quadratic  Regulator(LQR)  technique  is  used  for  control  computation.  Once  the 
synthetic  controller  v  is  evaluated  the  expression  for  the  actual  controller  of  the  plant  U  can  be 
evaluated  by  inverting  the  plant  dynamics.  Inverting  and  simplifying  the  above  approximate 
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plant  dynamics  (1.6)  the  expression  for  control  U  can  be  written  as 


U=  [, g{Xest)]~\v-f0iXest)-d{Xest)-Kr{Xest-Xa ))  (1.9) 

This  closed  form  expression  for  control  U  is  valid  provided  [g(Xesf)  ]  1  exists  for  all  values 
Xest.  In  particular  for  the  quad-rotor  plant  considered  in  this  paper  g(Xest )  is  a  product  of 
inverse  of  rotation  matrix  and  inertia  matrix  for  the  quad  frame.  Since  both  the  matrix  are 
invertible  at  all  times  except  when  roll  and  pitch  attitude  of  vehicle  i.e.  (f>,9  =  90°  the  closed 
form  expression  for  control  (1.9)  is  valid  for  all  other  domain  of  operation. 

1.1.2  Neural  Network  Synthesis  and  Weight  Update  Rule 

The  single  layer  NN  is  designed  to  capture  the  un-modeled  dynamics  of  the  plant.  It  is  desired 
that  virtual  (approximate)  plant  model  should  track  the  actual  plant.  Error  term  between 
actual  (X)  and  virtual  states  vector  ( Xa )  can  be  defined  as 

E  =  X-Xa  (1.10) 

Since  the  actual  state  vector  X  is  not  available  instead  a  EKF  estimate  is  used,  note  that  here 
we  introduce  a  small  abuse  of  notation  in  this  section,  state  vector  X  denotes  the  EKF  estimate 
of  the  state  vector.  The  error  dynamics  can  be  obtained  by  differentiating  the  above  equation 
with  time  and  substituting  for  X  and  Xa  from  (1.3)  and  (1.5) 

E  =  X-Xa 

E  =  d{X)-d{X)-KTE  (1.11) 

Single  layer  NN  with  RBF  basis  functions  is  chosen  to  approximates  the  un-modeled  dynamics 
dj  (X)  in  the  i th  channel. 


di(X)  =  W^dX),  Wi  e  9ftp  (1.12) 

where,  W)  are  the  actual  weights  and  ®,(X)  are  the  basis  function.  Lets  consider  there  exists 
an  ideal  approximator  for  the  unknown  function  which  approximates  d,  (X)  with  an  ideal 
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approximation  error  e;  for  the  chosen  basis  O,  (X). 


di(X)  =  WiT^i{X)  +  ei  (1.13) 

The  weights  Wj  are  the  ideal  weights  which  are  unknown.  Channel  wise  error  dynamics  can 
be  written  as 


et  =  Wi 7  ® (X)  +  e(-  -  Wi  7  0(-  (X)  -  kTi  et  (1.14) 

The  error  in  weights  of  the  i rh  approximating  network  is  defined  as 

Wi  =  Wi-Wi  (1.15) 

Wj  =  -Wi,  W{  =  constant  (1.16) 

Aim  is  that  weights  of  the  approximating  networks  Wj  should  approach  the  ideal  weights  Wj 
asymptotically,  i.e., 


Wi  —  0  as  t  —  oo 

The  un-modeled  information  is  stored  in  terms  of  the  weights  of  the  approximating  networks. 
Lyapunov  based  approach  for  deriving  weight  update  rule  is  discussed  in  the  next  subsection 
for  updating  Wj  online. 


Weight  Update  Rule 


The  candidate  Lyapunov  function  ensures  the  asymptotic  stability  of  the  following  variables, 
channel  wise  error  e, ,  error  in  network  weights  Wj  and  error  in  gradient  of  disturbance  term 
(for  directional  learning) 

ddj{X)  ddjUQ 
dX  dX 

A  positive  definite  Lyapunov  function  candidate  function  is  selected  as  follows 


=  WjT 


d®j 

dX 


(1.17) 


„  WjTWj 
Vi{ei}Wi)  =  Pi^-  +  — - - 

2  2  jj 


+  WjT 


d<Pj 

©/ 

dX 

2 

dX 

Wj 


(1.18) 
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where,  /3, ,  y,  and  0,-  are  positive  definite  quantities.  Taking  time  derivative  of  Lyapunov 
function  and  substituting  the  expression  for  error  dynamics  e,  from  (1.14)  and  simplifying 
the  expression  for  V,  can  be  written  as  follows 


V,-  = 


Pi e,-  Wjq>i  (X)  +  Pi ejel  -  p,  kTi  e] 


WjWi 


fid>; 


fid>; 


l 

Ti 

-Wi 

<5X 

6 

[  fix 

~  rr  (l 

+  WT  — 
dt 

d®i  1 
dX 

e 

fid),- 

fix 

T 

w 

w 


(1.19) 


Collecting  the  coefficients  of  W  and  equating  it  to  zero  the  following  weight  update  rule  in 
continuous  time  is  obtained 


P^i 


‘Ip 
—  + 
Ti 


d®j 

dXi 


0, 


<30; 


dXi 


r\_1  l 

<p{Xi)  + 


Ti 


a®. 


dXi 


0; 


dd>j 

dX< 


rx-i 

t 


X 


d 

dt 


a®; 

,  dX 


0 


30;  1 

— k 

dX 


T 

1.29) 


The  weight  update  rule  (1.20)  involves  the  term  W  on  RHS  of  the  expression,  i.e  to  propagate 
the  weight  matrix  it  is  required  to  have  the  information  of  term  Wjt  =  -  Wj.  which  in  turn 

demands  information  on  the  ideal  weights  at  all  time  step  ‘  k’.  Since  the  ideal  weights  are 
not  available  direct  information  of  the  term  W  cannot  be  extracted.  Having  said  that  the 
information  of  the  difference  of  the  ideal  weights  and  actual  weight  can  be  extracted  from  the 
difference  of  the  actual  disturbance  term  d(X)  and  approximated  disturbance  term  d[X).  It  is 
quite  evident  that  the  actual  disturbance  term  d{X)  is  unavailable  for  computation,  hence  an 
estimate  of  the  actual  disturbance  term  obtained  through  the  Extended  Kalman  Filter  (EKF)  is 
used  instead  of  the  actual  disturbance  term  d  (X)  in  the  above  expression  for  the  estimate  of 
the  term  W  the  details  of  which  are  given  in  subsequent  section. 

The  term  W  can  be  obtained  from  definition  of  the  term  (prW  as  follows 


cpTW=(dest(X)-d(X)) 


(1.21) 


note  that  the  actual  disturbance  term  is  replaced  by  the  EKF  estimate  of  the  d(X),  this  ap¬ 
proximation  is  reasonable,  since  the  estimate  of  d{X)  can  be  represented  in  terms  of  ideal 
weights  as  dest{X)  =  W  r<h(Xest).  Multiplying  through  out  by  the  term  0  and  rewrite  the  above 
equation  (1.21)  as 

(PT cpw  =  (p[dest{X)  -  d(X))  (1.22) 
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simplifying  the  above  equation  the  expression  for  the  term  W  is  obtained  as  follows 

W  =  l(pidestiX)  -  d(X))  (1.23) 

The  term  <pr<p  is  positive  definite  scalar  function  hence  the  inverse  always  exists,  substituting 
for  W  in  (1.20)  the  revised  weight  update  rule  is  obtained  as  follows. 


dO,- 

dO,- 

n 

-1 

(Ip 

dO,; 

dO,- 

r 

W  =  ptet\ 

— + 

dXi  \ 

dXi  \ 

J  (piXi)  +  | 

— + 

[dXi] 

[dXi\ 

[rt 

y 

d  /dO,t  ld®i\T  T  ,  ? 

d7(dx)0(dx)  x^)-V(^t(X)-d(X))  (i.24) 


(l„ 

dO,- 

dO,: 

T) 

— 1 

Piet 

— + 
r/ 

[ax/1 

0« 

dXt  \ 

) 

(ptXi)  + 

(< pT<p)-l(p({Xest  -  f[X)  -  bU)  -  d(X)) 


dP,- 

dXi 


®i 


dQ,- 

dXi 


d 

d  t 

,dX) 

dp\T 

dx) 


X 


(1.25) 


The  left  over  terms  from  Lyapunov  derivative  V,  equation  are 

h)  —  fiieiei  ~  kjifii 


(1.26) 


Using  Vi  <  0  leads  to  a  condition 


I  et\  > 


\£j\ 
kr  i 


(1.27) 


Therefore,  if  the  network  weights  are  updated  based  of  the  rule  given  in  (1.25),  then  the 
identification  happens  as  long  as  absolute  error  is  greater  than  t— .  By  increasing  kTj  error 

Kji 

bound  can  be  theoretically  reduced. 


1 .2  Generalized  Model  Predictive  Static  Programming 

In  this  section,  the  detailed  description  of  GMPSP  algorithm  is  presented.  General  form  of 
non-linear  state  dynamics  and  output  equation  is  defined  as, 

X(t)  =  f{X(t),U{t ))  (1.28) 
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Figure  1.1:  Feedback  Linearization  Trajectory  control  Augmented  with  RBF-NN,  control  Algo¬ 
rithm 


Y{t)  =  h{X(t ))  (1.29) 

Where  X  e  5R”,  U  e  5Rra  and  Y  e  5RP  With  known  boundary  conditions  X{to)  =  Xq  and  Y^itf)  as 
desired  output. 

The  philosophy  of  GMPSP  technique  is  to  use  error  history  between  present  and  desired  state 
value  and  update  control  history.  A  guess  control  history  is  required  to  start  GMPSP  algorithm 
which  yields  inaccurate  state  solution.  Error  history  is  computed  with  the  help  of  desired 
solution  and  inaccurate  solution  of  states.  This  error  history  is  use  in  update  continuously 
control  history  until  the  convergence  criteria  is  met  i.e.,  Y*  ( tf )  . 

Next,  mathematical  formulation  of  the  GMPSP  algorithm  is  presented,  Let  the  error  between 
present  and  desired  output  at  the  final  time  tf  be  defined  as  follows, 

6nX(tf))  =  [Yitf)-Y*{tf)]  (1.30) 

By  multiplying  both  sides  of  (1.28)  by  a  weighting  matrix  W{t)  produces 

w{t)x  =  w{t)f(xit),  mm  (i.3i) 

where,  the  computation  of  the  matrix  Wit)  e  SRPX"  is  described  later  (1.36)  in  this  section. 
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By  integrating  both  sides  of  (1.31)  from  to  to  tf  leads  to, 


tf  tf 

J  [wmx){t)]dt  =  J  [ W(t)f(X(t),U(t))]dt 


(1.32) 


Next,  on  adding  the  quantity  Y ( X[tf ))  to  both  sides  of  (1.32)  and  by  using  algebraic  manipula¬ 
tion,  the  following  is  obtained, 


*7 

■/ 


Y{X(tf))  =  Y{X(tf))+  [ W{t)f{X(t),mt))]dt -  [ W(t)X{t)]dt 


V 

/■ 


(1.33) 


Using  the  integration  by  parts  in  the  last  term  of  the  right  hand  side  of  (1.33)  ,  the  following 
can  be  written: 


Y{X{tf))  =  Y{X{tf))  -  [W{tf)X{tf]\  +  [W{to)X{to) 

d 

+  j  [Wit)f{X(t),mt))  +  W(t)X(t)]dt  (1.34) 


By  taking  the  first  variation  of  the  both  sides  of  (1.34), 


SY{X[tf))  = 


' dYVOj )) 
,  dX(t) 


—  W{t)\  5X{t) 


t=tf 


rt  df[X(t),mt))  .  p 

+  \W[t)  J  — d- - +  W{t)  SXWdit) 

J  l  dX{t) 


+  [Vfi(fo)<5X(fo)]  + 


djmtwuitiv 

dU{t) 


8U[t)d{t) 


(1.35) 


Next,  it  is  desired  to  determine  the  variations  8Y ( X(tf ))  produced  by  the  given  8U(t)  .  The 
idea  is  to  choose  W(t)  in  such  a  way  that  vanishes  the  coefficients  of  5X(t)  in  the  above 
equation.  Weighting  matrix  W[t)  can  be  selected  as, 


W(t)  =  -W(t) 


dX{t ) 


(1.36) 
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W(tf)  = 


dY{Xitf )) 
dX{if ) 


(1.37) 


Initial  condition  is  known  Hence,  the  expression  SXito)  =  0  holds  true  .  Furthermore,  substi¬ 
tuting  (1.36)  and  (1.37)  into  (1.35)  produces, 


* f 

8Y(X{tf))  =  j  [Bs{t)8U{t)]dt 

to 


Where, 


Bs(t)  =  Wit) 


dfjXitlUm 

dUit) 


For  optimal  control  formulation,  cost  function  is  defined  as, 


J'  =  \\  mp(t)-8int))TR(t){up(t)-su{t))]dt 

to 


(1.38) 


(1.39) 


(1.40) 


where,  R{t)  >  0  is  the  non-singular  matrix.  The  cost  function  in  (1.40)  needs  to  be  minimized, 
subjected  to  the  constraint  in  (1.38).  With  the  application  of  optimization  theory  Bryson  and 
Ho  [1975],  static  constrained  optimization  problem  can  be  formulated  by  using  Equations 
(1.38)  and  (1.40).  The  augmented  cost  function  is  given  as, 

f/ 

j=U  ttupw-5uwTR^upw-5u^dt 

to 

lf 

+AT[8Y(tf)-J  [Bs(t)8U(t)]dt]  (1.41) 

to 

Where,  A  is  the  Lagrange  multiplier.  First  variation  of  (1.41)  leads  to, 

8j=- J  mmup(t)-smt))  +  (Bs{t))TA8i8uim  (1.42) 

to 

From  (1.42)  it  is  clear  that  a  minimum  of  /  occurs  if  the  following  expression  holds  true: 


SUit)  =  iRit))~1iBsit))TA  +  Upit) 


(1.43) 


Distribution  Code  A:  Approved  for  public  release,  distribution  is  unlimited 


The  closed  form  solution  of  control  history  U{  t)  is  obtained  using  GMPSP  algorithm.  The  basic 
mechanism  behind  GMPSP  technique  is  to  convert  the  dynamic  optimal  control  problem 
into  a  constrained  static  optimization  problem  and  update  the  control  history  using  small 
dimension  weighting  matrix  update  (1.36).  This  leads  to  fast  computation  of  control  history. 
Due  to  fast  convergence  property  GMPSP  algorithm  can  be  applied  for  online  applications. 
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Autonomous  Soft  Landing  of  Quadro- 
tor 


Owing  to  highly  agile  nature  of  the  quadrotor  vehicle  and  inherent  existence  of  exogenous 
disturbing  forces,  safe  landing  of  the  quadrotor  with  safe  approach  velocities  is  always  a  tricky 
and  critical  issue  for  safety  of  the  vehicle  and  surroundings.  Landing  scenario  becomes  still 
more  critical  when  safe  landing  is  to  be  achieved  on  the  specified  location  within  a  confined 
space  with  obstacles.  For  a  non  qualified  pilot,  safe  landing  of  the  vehicle  is  a  challenging 
task  to  accomplish.  External  disturbances  like  cross  wind,  blade  flapping,  ground  effect 
on  thrust  variation  due  to  downwash  vehicle  dynamics  change  according  to  possible  flight 
scenarios.  Apart  from  the  external  factors  measurement  noise  also  hinders  safe  decent  through 
error  in  position  and  attitude  measurement  of  the  vehicle.  Vision-based  methods  for  aerial 
vehicles  to  detect  navigation  targets  have  been  used  recently  and  reported  in  literature  to 
identify  the  targets  for  path  following  and  landing.  Brockers  et  al.  Brockers  et  al.  [2011]  have 
presented  Autonomous  landing  and  ingress  of  micro -air- vehicles  in  urban  environments 
using  monocular  vision  processing.  Templeton  et  al.  Templeton  et  al.  [2007]  demonstrated  a 
model-predictive  flight  control  system  using  a  monocular  camera.  3D  feautures  are  tracked  in 
earth  fixed  inertial  frame  to  estimate  elevation  and  appearance  of  the  structure,  their  approach 
assumes  the  availability  of  GPS  information  for  computation.  Recently  Blosch  et  al.  Blosch 
et  al.  [2010]  highlight  in  their  work  the  use  of  a  monocular  vision-based  SLAM  algorithm  for 
accurate  determination  of  the  pose  of  a  UAV  in  GPS  denied  environments.  This  work,  while 
not  exploring  techniques  landing  site  selection  per  se,  dwells  more  towards  the  nonlinear 
controller  and  online  adaptive  disturbance  estimation  techniques  augmenting  the  capability 

15 


Distribution  Code  A:  Approved  for  public  release,  distribution  is  unlimited 


2.1  Quadrotor-  Plant  Dynamics 


16 


of  operating  UAV  systems  in  all  types  of  environments  and  concentrate  on  the  autonomous 
landing  capability  of  the  vehicle.  For  simulations,  the  problem  statement  of  landing  at  location 
X  =  5m,  Y  =  5m  and  Z  =  0.05 m  assumed  to  be  identified  by  available  ranging  techniques  is 
used.  The  vehicle  is  assumed  to  be  hovering  above  ground  level  at  Z  =  10m  a  error  in  initial 
condition  in  X  and  Y  direction  is  assumed,  the  vehicle  is  off  by  lm  in  each  direction  that 
is  X  =  4m,  Y  =  4m.  From  this  initial  position  of  the  quadrotor  the  controller  is  designed  to 
achieve  auto-landing  at  predetermined  location,  under  the  influence  of  external  disturbance 
of  rotor  blade  flapping,  propeller  axis  misalignment.  The  control  architecture  and  simualtion 
results  are  presented  in  the  further  section. 


2.1  Quadrotor-  Plant  Dynamics 

This  section  provides  the  details  of  kinematics  and  dynamic  equation  of  motion  of  the  rigid 
body  quadrotor.  The  notation  and  the  co-ordinate  frames  used  are  typical  in  aeronautics 
literature,  for  further  details  one  can  refer  Beard  [2008].  Euler  angles  in  sequence  of  1-2-3 
(Roll  0,  Pitch  6  and  Yaw  \ fr)  is  used  for  deriving  the  transformation  matrix  for  translation  and 
angular  velocity  for  the  vehicle.  This  chosen  frame  for  transformation  has  advantage  that  it 
renders  the  outer-loop  control  i.e.  position  control  independent  of  yaw  attitude  of  the  vehicle. 
Quadrotor  Dynamics  equation  of  motion  is  expressed  as  follows 


f  = 

(p  = 


+ 


-G+—R,F 

M 

5-1  7-1 


-R-/F1  [Rr(p  x  JRr(p) 


R 


-l 


3Rr  •  dRr  . 


Rr  J  T  +  dexternai{X) 


(2.1) 


(2.2) 


where  rotation  matrices  for  transforming  translation  and  angular  position  of  the  Quadrotor 
from  body  frame  to  inertial  frame  are  as  follows 


Cq  Cy/ 

—  C()Sy/ 

Se 

S(p  Sq  Cy/  +  C0  Sy 

—  Sff)  Sq  Sy/  +  C0  Cy/ 

-S^Cg 

—  C(f)  Sq  Cy/  +  S(f)  Sy/ 

C(ff  Sq  Sy/  +  S(f)  Cy/ 

C(pCg 
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Rr  = 


Cy/  Cg  Sy/  0 
—  Sy,  Cg  Cy/  0 

Sg  0  1 


and  vectors  cf  =  [jc,  y,z] T  and  cp  =  [<p,d,y/]T  are  position  and  angle  vector  of  the  quadrotor 
in  inertial  frame  of  reference.  V  =  [u,v,w]T  and  O  =  [p,  q,  r] T  are  translation  velocity  and 
rotation  velocity  in  body  fixed  reference  frame  and  G  =  [0, 0,  g]  T  denotes  gravity  vector.  Inter¬ 
ested  readers  can  refer  Beard  [2008],  Lee  et  al.  [2011]  and  references  there  in  for  the  detailed 
derivation  of  the  equation  of  motion  of  the  Quadrotor  dynamics. 


2.1.1  Quadrotor  Blade  Flapping  Phenomenon 

A  quadrotor  UAV  in  linear  translational  flight  experiences  a  aerodynamic  effect  known  as 
“Blade  Flapping  Phenomenon”.  This  aerodynamic  phenomenon  causes  disturbing  moments 
acting  about  the  c.g  of  the  vehicle.  The  disturbing  moments  are  caused  as  the  result  of  the 
unequal  lifting  forces  experienced  by  the  rotor  blade  due  the  variation  in  the  free  stream 
velocity  experienced  by  the  leading  edge  and  trailing  edge  of  the  rotor  blade.  In  steady  state 
this  effect  causes  a  steady  bias  in  the  rotor  blade  angle  causing  tilt  in  thrust  vector  and  a 
component  of  the  trust  vector  parallel  to  Roll-pitch  plane  causes  the  moments  acting  at  the 
c.g  of  the  vehicle.  The  steady  bias  in  the  rotor  blade  tilt  due  blade  flapping  phenomenon 
also  causes  a  steady  moment  due  structural  deflection  of  the  blade,  this  effect  is  observed 
in  the  more  stiff  propellers,  the  reaction  moment  at  the  rotor  hub  is  product  of  the  stiffness 
coefficient  of  the  rotor  blade  and  deflection  caused  due  to  flapping. 

The  tilt  of  the  rotor  plane  through  angle  ^  in  pitch  and  £,g  in  roll  directions,  generates  a  thrust 
vector  component  in  longitudinal  and  lateral  directions  causing  disturbing  moments  about 
pitch  and  roll  axis, 


M pitch  —  Thsin^y)  +  Up  t  y> 
Mr  oil  =  Thsin{g  +  Kptg 


(2.3) 

(2.4) 
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where  T  is  total  thrust  (N),  h  is  the  perpendicular  distance  between  rotor  plane  and  roll-pitch 
plane  passing  through  c.g  of  the  vehicle  and  £,,p,£,o  are  the  rotor  plane  deflection  angles. 

The  rotor  plane  deflection  angle  is  function  of  the  translation  velocities  of  the  vehicle  w.r.to 
the  free  stream  of  air  in  inertial  frame  and  roll  and  pitch  attitude  angular  rates  of  the  quadrotor 
vehicle  and  can  be  evaluated  as  follows 


t*  = 

C\X  +  C20 

(2.5) 

te  = 

Ciy+ C2O 

(2.6) 

where  Ci  and  C2  are  aerodynamic  constants  depending  on  the  rotor  blade  properties.  Full 
analysis  of  the  blade  flapping  phenomenon  is  beyond  the  scope  of  this  paper,  interested  reader 
can  refer  Huang  et  al.,  Hoffmann  et  al.  and  references  there  in  for  further  details. 

2.2  Neuro-Adaptive  Dynamic  Inversion  for  Quadrotor  Control 

The  baseline  controller  is  evaluated  using  dynamics  inversion  method  Lee  et  al.  [2011].  The 
Guidance /Outer  loop  dynamics  concerns  the  position  control  of  the  vehicle.  The  linear 
feedback  representation  of  outer-loop  system  dynamics  can  be  written  as 

£  =  w 

Linear  control  technique  such  as  linear  quadratic  regulator  (LQR)  is  used  to  evaluate  the 
control  w.  Control  w  is  of  form  w  -  -Kposition(£  -  £,d).  Now  using  synthetic  controller  w  the 
total  thrust  vector  is  evaluated  as  follows 

1 

w  =  -G+—RtF  (2.7) 

M 

denoting  U  =  RtF  from  above  equation  U  can  be  expressed  as  U  =  Mw  +  MG  using  the 
definition  U  =  RtF  the  expression  for  total  thrust  F  is  solved  as  follows 

F=Jiu[~Gufcu[) 
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and  desired  roll  and  pitch  attitude  is  evaluated  as  follows 

r 

(pd  =  -sgn{U2)cos~l 

k 

and 

r 

0d  =  -sgn[Ui)cos~l 

k 

Desired  yaw  angle  is  considered  to  be  y/d  =  0.  The  desired  attitude  angles  forms  input  to  the 
innerloop  controller  of  the  quadrotor. 

Similarly  the  inner  loop  dynamics  can  be  written  in  linear  form  for  the  feedback  linearization 
controller  evaluation  as  follows. 


U23 


u2  +  ui 


T  =  JRr{v  +  R-rlrl{Rr(jpxJRrq))) 


+  JR, 


r: 


3Rr 

{38 


■  3Rr  . 
■6+—y/ 
dyr 


(p  +  d{X) 


(2.8) 


2.3  Extended  Kalman  Filter(EKF)  Estimate  of  the  Quadrotor  State 
Vector 

The  quadrotor  state  vector  is  available  through  the  sensor  measurement.  On-board  sensors 
like  three  axis  accelerometers  and  rate  gyro  provides  the  body  rate  information  with  respect  to 
the  body  frame  and  GPS  measures  the  position,  ground  speed  and  course  angle  with  respect  to 
inertial  frame.  Rate  gyros  and  accelerometers  output  body  rate  information  and  is  assumed  to 
be  corrupted  by  zero  mean  measurement  noise  and  process  noise  arising  due  to  blade  flapping 
phenomenon,  hence  a  EKF  estimate  needs  to  be  propagated  for  information  on  the  attitude 
state  vectors.  Where  as  the  position  and  velocity  of  the  quadrotor  in  inertial  frame  is  assumed 
to  be  faithfully  available  through  GPS  measurement.  Hence  the  EKF  is  only  implemented 
for  the  inner-loop  dynamics,  there-by  with  directly  available  outerloop  state  information 
the  EKF  estimate  of  inner  loop  state  vector  forms  the  complete  state  information  of  the 
quadrotor  vehicle.  Fig  2. 1  provides  details  to  EKF  philosophy  implementation  for  quadrotor 
state  estimation,  for  further  details  interested  readers  can  refer  Crassidis  and  Junkins  [2011]. 
The  estimate  of  the  uncertainty,  derived  from  EKF  estimate  of  the  innerloop  state  vectors  for 
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Figure  2.1:  Extended  Kalman  Filter  Implementation  for  Partial  State  Estimation  (Inner  Loop 
Attitude  State  Estimates) 


weight  update  rule  (1.20)  is  evaluated  as  follows 

d(Xest)  =  Xest  -  fo(X)  -  g{X)U  (2.9) 


2.4  Numerical  Simulation  for  Quadrotor  Landing 

The  estimate  of  the  uncertainty,  derived  from  EKF  estimate  of  the  innerloop  state  vectors  for 
weight  update  rule  (1.20)  is  evaluated  as  follows 


d(Xest)  =  Xest  -  MX)  -  g{X)  U 


(2.10) 


2.4.1  Simulation  Results 

This  section  presents  simulation  results  for  quadrotor  Landing  control.  The  control  law  de¬ 
signed  treats  the  disturbing  moment  due  to  aerodynamic  effect  of  blade  flapping  phenomenon 
as  un-modeled  uncertainty  term  which  is  captured  online  through  RBF-NNs.  Outer  loop  con¬ 
trol  is  achieved  through  a  dynamic  inversion  control  and  inner  loop  attitude  control  is  achieved 
through  the  online  adapted  dynamic  inversion  control.  At  low  speeds  and  less  aggressive 
maneuver  the  aerodynamic  effects  due  quadrotor  blade  flapping  phenomenon  are  less  sever 
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Quad -Rotor  Trajectory  Tracking 


Figure  2.2:  Quad  Rotor  following  the  pre  Defined  Trajectory:  Reference  Trajectory,  Actual  and 
Approximate  trajectory  of  vehicle 


Quad-Rotor  Roll  Attitude 


Figure  2.3:  Reference/Desired,  Actual,  Measured  and  Kalman  filter  Estimate  of  die  Roll  Attitude 
of  the  Quad 
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Quad-Rotor  Pitch  Attitude 


Figure  2.4:  Reference /Desired,  Actual,  Measured  and  Kalman  filter  Estimate  of  the  Pitch 
Attitude  of  the  Quad 


Quad-Rotor  Yaw  Attitude 


Time  (sec) 


Figure  2.5:  Reference/Desired,  Achieved  Yaw  Attitude  of  the  Quad 
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and  it  is  seen  that  a  well  tuned  PID  controller  is  sufficient  for  satisfactory  command  tracking  In 

Error  in  Attitude  Angles 

0.1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 


-e- 


Figure  2.6:  Error  in  Attitude  Angle  of  the  Quadrotor 


Figure  2.7:  Quad  Rotor  following  the  pre  Defined  Trajectory:  Reference  Trajectory,  Actual  and 
Approximate  trajectory  of  vehicle 


translational  flights  it  is  observed  that  at  larger  speeds  the  blade  flapping  causes  a  considerable 
amount  of  disturbing  moment  at  the  c.g  of  the  vehicle.  At  higher  speeds  the  restoring  forces 
may  not  be  sufficient  and  vehicle  eventually  becomes  unstable  unless  the  control  architecture 
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Disturbance  Approximation 


Figure  2.8:  Disturbance  in  Channel  1,  Actual,  Approximated,  Kalman  Filter  Estimated 


implemented  takes  into  account  the  blade  flapping  phenomenon.  The  performance  of  the 
proposed  controller  in  auto-landing  of  the  vehicle  from  altitude  of  10m  is  shown  in  the  Fig 
2.2.  The  desired  location  of  the  landing  is  X  =  5m,  Y  =  5m  and  Z  =  0.05m.  The  vehicle  initial 
condition  at  point  of  start  of  decent  for  landing  is  considered  t  be  at  X  =  4m,  Y  =  4m  and 
Z  =  10m  above  ground  level.  Quadrotor  roll  attitude  against  the  desired  commanded  attitude 
is  shown  in  the  Fig  2.3.  Figure  2.3  also  gives  the  details  of  the  measured  attitude  angle  and  EKF 
estimate  of  the  roll  attitude  of  the  vehicle. 


Disturbance  Approximation 


Figure  2.9:  Disturbance  in  Channel  3,  Actual,  Approximated,  Kalman  Filter  Estimated 


Similarly  Fig  2. 4, 2.5  gives  the  details  of  the  attitude  control  of  the  quadrotor  in  pitch  and 
yaw  axis  respectively.  It  is  observed  that  despite  the  influence  of  the  disturbing  aerodynamic 
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Disturbance  Approximation 
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Figure  2.10:  Disturbance  in  Channel  3,  Actual,  Approximated,  Kalman  Filter  Estimated 


Propellor  Thrust  Magnitude 


Figure  2.11:  Control, Thrust  Produced  by  four  propellers 


moments  the  perfect  trajectory  tracking  is  achieved.  The  vehicle  attitude  follows  the  desired 
attitude  angles  with  high  accuracy,  the  errors  in  the  attitude  angles  is  presented  in  Fig  2.6.  The 
plot  of  state  history  for  the  commanded  trajectory  following  in  X,  Y  and  Z  in  inertial  frame  of 
reference,  are  provided  in  the  Fig  2.2. 

The  performance  of  the  three  RBF-NN  in  approximating  the  disturbance  terms  in  each 
channel  is  given  by  the  Fig  2.8, 2.9  and  2.10.  Figure  2.8, 2.9  and  2.10  provides  the  details  of 
the  actual  disturbing  moments,  EKF  estimated  disturbance  and  neural  network  approximate 
of  the  uncertainty  arising  due  to  blade  flapping  phenomenon.  It  can  be  observed  that  EKF 
estimation  of  d(X)  through  estimation  of  acceleration  term  Xest  leads  to  amplification  of  the 
effect  of  the  sensor  noise  into  the  uncertainty  estimate.  Authors  wishes  to  highlight  the  notable 
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Figure  2.12:  Weight  history,  W1 


Figure  2.13:  Weight  history,  W2 


advantage  of  weight  update  rule  presented  in  this  paper,  using  the  noisy  EKF  estimate  of  the 
uncertainty  into  the  weight  update  rule  (1.25)  a  precise  approximation  is  achieved  through  the 
neural  network  approximation  for  unknown  system  dynamics.  Figure  2.11  provides  the  details 
of  the  final  control  effort  needed  in  terms  of  the  propeller  thrust  to  achieve  the  commanded 
trajectory  tracking.  The  control  effort  F\ ,  F2,  F3,  F4  corresponds  to  thrust  generated  by  front, 
right,  back  and  left  propellers  respectively.  Figure  2.12,  2.13  provides  time  history  plot  of  the 
weight  propagation  for  online  approximation  of  the  uncertainty.  The  value  of  weight  W3  are 
negligible  since  the  disturbance  in  yaw  axis  is  zero  hence  the  plot  is  not  provided. 
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Figure  2.14:  Weight  history,  W3 
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Autonomous  Soft  Landing  of  Lunar 
Module 


The  powered  descent  phase  starts  from  perilune  as  the  nearest  point  from  the  Moon  surface. 
During  the  powered  descent  phase  the  reverse  thrust  action  provides  a  braking  mechanism  to 
reduce  the  high  orbital  velocity  (about  1700m/.?)  of  the  spacecraft.  At  the  end  of  the  powered 
descent  phase  the  spacecraft  is  positioned  very  close  to  the  lunar  surface  about  200  m  altitude 
with  a  reduced  velocity  in  the  range  of  (0  -  20)  m/s. 


2 


T 


X 


Figure  3.1:  Landing  site  based  frame  of  reference 


3.1  Problem  Statement 

Due  to  the  uneven  terrain  of  lunar  surface,  the  visual  sensors  are  used  to  detect  a  suitable 
landing  site.  It  might  be  necessary  to  shift  the  landing  site  from  the  nominal  location  to  a 
nearby  smooth  ground.  The  operation  of  switching  the  landing  site  is  known  as  re-targeting. 
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Hence  it  is  necessary  to  obtained  an  on  board  guidance  law  that  can  ensure  the  soft  landing 
requirement  while  accounts  for  the  correction  of  nominal  landing  site.  An  analytical  guidance 
law  for  fuel  optimal  planetary  landing  was  studied  in  D’Souza  [1997] .  A  close  form  expression 
of  optimal  commanded  acceleration  has  been  obtained  for  the  spacecraft  equation  of  motion 
in  terms  of  current  position  and  velocities  and  also  a  function  of  the  time-to-go  parameter.  The 
guidance  law  is  popular  as  constrained  terminal  velocity  guidance  (CTVG) .  A  relations  between 
the  CTVG  with  the  classical  proportional  navigation  (PN),  zero  effort  velocity  miss  (ZEVM) 
and  zero  effort  miss(ZEM)  feedback  guidance  laws  are  presented  in  Guo  et  al.  [2011].  The 
CTVG  guidance  ensures  the  soft  landing  but  fails  to  demonstrate  the  vertical  orientation  of  the 
spacecraft  during  the  terminal  descent.  A  modified  CTVG  guidance  considers  the  derivative 
of  the  acceleration  command  as  a  virtual  input  vector.  The  terminal  attitude-constrain  is 
incorporated  in  the  modified  CTVG  law,  is  presented  in  Zhao  et  al.  [2014].  The  computation  of 
the  modified  time  to  go  presented  in  Zhao  et  al.  [2014]  is  a  very  complex  recursive  relation. 
Authors  have  investigated  the  numerical  method  for  the  optimal  control  design  for  the  lunar 
landing.  Zhang  Jianhui  Jianhui  et  al.  [2011]  uses  the  Chebyshev  pseudo-spectral  method  for 
optimal  control  design  using  two  dimensional  planar  motion  of  lunar  module.  Time  scaling 
transformation  followed  by  Control  vector  parameterization  is  presented  in  Gao  et  al.  [2013] , 
Liu  et  al.  [2008]  as  trajectory  optimization  method  for  the  planner  dynamics  of  lunar  module. 
Although  the  various  numerical  optimal  control  methods  are  studied  for  optimal  guidance 
of  lunar  modules, few  are  suitable  for  onboard  implementation  due  to  poor  computational 
efficiency  with  low  convergence  rate. 

G-MPSP  is  a  numerical  optimal  control  algorithm  deals  with  general  non-linear  systems.  The 
algorithm  is  based  on  the  philosophies  of  the  approximate  dynamic  programming  and  nonlin¬ 
ear  model  predictive  control.  The  formulation  of  G-MPSP  converts  the  dynamic  optimization 
problem  (in  HJB  framework)  into  a  static  programming  problem  results  a  close  form  solution 
for  the  control  update  history.  The  close  form  control  update  law  makes  the  G-MPSP  com¬ 
putationally  efficient.  The  fast  convergence  rate  makes  it  compatible  for  on-board  optimal 
guidance  design.  In  this  report  G-MPSP  guidance  is  formulated  for  the  accurate  soft  landing  of 
a  lunar  module.  From  the  mission  point  of  view,  the  module  starting  from  the  altitude  about 
200  m  requires  to  touchdown  with  near  zero  velocity  about  2  m  altitude  from  the  designated 
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landing  site  (determined  by  the  vision  sensor)  on  lunar  surface.  For  a  spacecraft  with  variable 
thrust  engine, G-MPSP  generates  the  acceleration  vector  (both  magnitude  and  direction)  that 
minimizes  the  fuel  requirement  and  ensures  the  accurate  soft  landing  of  the  module. 


3.2  Spacecraft  Dynamics  for  Lunar  Mission 

The  spacecraft  equation  of  motion  derived  in  the  landing  site  based  inertial  frame  of  reference 
is  described  as, 


x=  u 


y=v 
z  =  w 


ui  T 

u  = - — -  H - cos  pcoscp 
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m  =  -- 
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(3.1) 


(3.2) 


where,  (x,  y,  z)  topple  represents  the  instantaneous  position  of  the  spacecraft  in  the  inertial 
frame  of  reference,  u,  v  and  w  denotes  the  velocity  component  of  the  module.  The  thrust 
engine  is  mounted  on  the  spacecraft  body,  Tim  describes  the  acceleration  vector  generated 
using  the  thrust  engine.  The  direction  of  the  thrust  vector  in  3-D  space  is  represented  by  the 
angle  /)  and  cp.  During  the  terminal  descent  the  objective  is  to  guide  the  spacecraft  from  an 
altitude  of  about  200 m  to  the  designated  location  on  the  lunar  surface.  The  lunar  gravity 
parameter  p  is  constant  for  the  domain  of  operation.  The  initial  position  and  the  velocity  of 
the  spacecraft  is  considered  as,  (xo,yo,-Zo,  uo,  vq,  w)-  The  specific  impulse  of  the  thrust  engine  is 
represented  by  Isp  and  go  is  the  Earth  gravitational  acceleration. 
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3.3  Fuel  Optimal  Guidance 

The  guidance  objective  is  to  obtain  the  thrust  vector  (both  magnitude  and  direction)  T,  p  and 
cp  such  that  the  consumption  of  the  fuel  mass  is  minimized,  so  the  cost  function  J\  is  defined 
as, 


Ja  ~ 


mo 

mf 


(3.3) 


The  initial  mass  of  the  spacecraft  m0  is  known  and  constant.  Hence  the  minimization  of 
Ja  is  equivalent  of  maximization  of  the  spacecraft  landed  mass.  The  minimization  of  Ja  is 
mathematically  equivalent  to  minimization  of  the  following , 


Jb  =  In  Ja  =  In 


'mo' 
'  mf  ‘ 


(3.4) 


Using  the  mass  dynamics  given  in  (3.1)  the  above  expression  can  be  represented  as, 


In 


' mo] _ _  r 
,mf)  J 


d{\nm(t)) 

dt 


dt 


Jb  =  ~ 


dt 


(3.5) 


By  using  (3.1) ,  (3.5)  leads  to, 


Now  the  acceleration  components  can  be  written  as, 


(3.6) 
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ax  =  — cos  Bcos(p 
m 


fly  =  — cosdsincp 
m 

T 

az  =  — sinp 
m 


(3.7) 


by  substituting  ax,  ay  and  az  in  (3.6)  one  can  obtain, 

ff 

Jb  =  J  \J cix  +  a2  +  a2zdt  (3.8) 

o 

here,  Jb  is  a  convex  function  with  respect  to  ax,  ay  and  az.  Minimization  of  //,  is  equivalent  to 
minimization  of 

J  =  J  ^a2  +  a2  +  a2^dt  (3.9) 

o 

Therefore, maximization  of  the  spacecraft  mass  of  the  lunar  module  during  landing  is  ensured 
by  minimizing  (3.9),  which  results  into  minimization  of  fuel  consumption  during  terminal 
descent  phase. 

3.4  Results  and  Discussion 

The  numerical  values,  considered  for  the  simulation  purpose  are  given  as  below.  The  nominal 
values  of  the  initial  state  variables  are  considered  as,  xo  =  0  m,  yo  =  0  m,  zq  =  200  m  represents 
the  initial  position  of  the  module  with  respect  to  the  landing  site  frame  of  reference,  uq  = 
5 ml  sec,  uq  =  -5ml  sec  and  wq  =  5 ml  sec  describes  the  initial  velocity.  The  initial  mass  of  the 
module  at  perilune  M0  =  523.5 kg.  The  objective  is  to  generate  the  thrust  vector  such  that  the 
module  can  land  safely  on  a  designated  landing  site  determined  using  the  visual  sensors.  At 
the  end  of  the  terminal  descent  the  state  variables  are  required  to  be  obtained  as, 
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Figure  3.2:  Spacecraft  trajectory 


Figure  3.3:  Position  co-ordinate  of  the  lunar  module 


Xf  =  100  m,  yf  =  20  m  and  Zf  =  0  m  (the  co-ordinate  of  the  detected  landing  site  in  the  same  in¬ 
ertial  frame  of  reference).  Uf  =  0ml  sec,  Vf  =  0ml  sec  and  Wf  =  0m /  sec  denotes  the  terminal 
touch  down  velocity.  The  lunar  gravity  parameter  p  =  4.90278  x  103  km3  /  sec2 ,  Isp  =  315  sec. 
Control  weight  matrix  Rjdt)  is  selected  as  a  time  varying  matrix  to  ensure  the  vertical  landing 
of  the  module.  The  G-MPSP  based  terminal  descent  guidance  generates  the  optimal  com¬ 
manded  thrust  for  the  lunar  module.  Fig.3.2  and  Fig.3.3  shows  the  spacecraft  (point  mass) 
trajectory  is  starting  from  the  specified  initial  position  and  follows  a  smooth  curve  and  ap¬ 
proaches  towards  the  detected  landing  site.  At  the  final  flight  time  the  distance  of  the  lunar 
module  from  the  ground  surface  is  about  z  —  Zf  =  2.  At  the  altitude  of  2m,  the  thrust  engine  of 
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Time  (sec) 

Figure  3.4:  Velocity  components  of  the  lunar  module 


Figure  3.5:  Velocity  vector  of  the  lunar  module 


the  spacecraft  remains  switch  off.  Hence  from  2m  altitude  the  spacecraft  will  be  allowed  to 
free  fall  until  it  touches  the  ground.  The  2m  margin  is  necessary  while  the  spacecraft  touches 
the  lunar  surface,  as  the  thrust  engine  plume  can  make  a  dusty  environment  over  the  solar 
panels  mounted  on  the  lunar  module.  To  ensure  the  soft  landing  as  the  time  approaches,  the 
spacecraft  velocity  approaches  zero  as  shown  in  Fig.3.4  and  Fig.3.5.  Starting  from  a  near  zero 
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Figure  3.6:  Thrust  direction  beta  profile 


Figure  3.7:  Thrust  magnitude  profile 


initial  guess  control  history  the  guidance  algorithm  takes  four  to  five  iterations  to  converge. 
The  optimal  commanded  thrust  vector  generated  using  G-MPSP  guidance  for  the  variable 
thrust  engine  is  shown  in  the  Fig.3.7  and  Fig.3.6. 
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Time  (sec) 


Figure  3.8:  Variation  spacecraft  total  mass 


The  variation  of  the  spacecraft  mass  is  shown  in  Fig.3.8,  the  fuel  requirement  during  the 
terminal  descent  phase  is  about  14A;g. 
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CONCLUSION 


Tow  problems  are  investigated  in  this  report.  First,  the  autonomous  guidance  of  soft  landing 
for  an  UAV  is  presented.  In  presence  of  atmosphere  with  the  utility  of  the  aerodynamic 
advantage  the  quadrotor  landing  on  the  flat  earth  surface  is  demonstrated  with  a  neuro- 
adaptive  tracking  control  approach  for  autonomous  landing.  The  approach  leads  to  fast 
learning  of  the  disturbance  function  in  the  system  dynamics  and  with  lesser  transients.  The 
presented  adaptive  control  was  experimented  with  quad-rotor  autonomous  landing  problem 
in  presence  of  aerodynamics  uncertainty.  It  has  been  shown  that  the  proposed  approach  could 
satisfactorily  capture  the  uncertainty.  The  inner-loop  neuro-adaptive  dynamic  inversion 
controller  forms  a  very  robust  controller  and  is  demonstrated  to  be  tracking  the  desired 
attitude  angles  with  high  accuracy  in  presence  of  unmodeled  dynamics  and  measurement 
sensor  noise.  As  a  result  the  outer-loop  controller  is  able  to  perfectly  track  the  commanded 
trajectory. 

The  second  problem  deals  with  the  autonomous  soft  landing  of  a  lunar  module.  In  absence 
of  atmosphere,  where  the  aerodynamic  advantages  are  not  present  (which  is  true  for  moon), 
the  terminal  descent  of  the  autonomous  guidance  for  a  spacecraft  on  the  flat  moon  surface  is 
demonstrated  using  the  recently  developed  G-MPSP  algorithm.  The  proposed  G-MPSP  guid¬ 
ance  ensures  the  objective  of  the  spacecraft  soft  landing  along  with  the  minimum  propellant 
consumption  with  high  terminal  accuracy.  With  the  proper  choice  of  the  control  weight  matrix 
the  terminal  vertical  orientation  of  the  spacecraft  is  obtained  in  the  soft  constraint  framework. 
The  effectiveness  of  the  proposed  algorithm  has  been  demonstrated  using  simulation  results. 
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